\documentclass{article}
\usepackage{amsmath}
\usepackage{graphicx}

\begin{document}
\section{Shan-Chen Multiple component model}
\label{chapterShanchen} A multi-component system consist of separate
chemical components such as oil and water.Due to their economic
importance, such systems has been studied extensively.\\

Adding a second chemical component without long range interactions
enable us to simulate completely miscible fluid, while completely
immiscible fluid simulation is possible where long range repulsive
interactions are considered in the model.\\

We add a second component by introducing a new index to the relevant
distribution functions: $f^{\sigma}$ represents the distributiong
functions of component $\sigma$. Similar definitions apply to
velocities, densities. To conserve the momentum locally in the
collision step, a composite macroscopic velocity is defined as:

\begin{equation}
u'=\frac{\sum\limits_{\sigma} \frac{1}{\tau_{\sigma}} \sum\limits_a
f_a^{\sigma}
e_a}{\sum\limits_{\sigma}\frac{1}{\tau_{\sigma}}\rho_{\sigma}}
\end{equation}

Uncoupled velocities and densities of individual component $\sigma$
are defined as usual:


\begin{equation}
\rho_a = \sum\limits_{a=0}^8 f_a^{\sigma}
\end{equation}

\begin{equation}
u_{\sigma}=\frac{1}{\rho_{\sigma}} \sum\limits_{a=0}^8 f_a^{\sigma}
e_a
\end{equation}

The force applied on component $\sigma$ is
\begin{equation}
F_{\sigma}(x)=-G \psi_{\sigma}(x,t)\sum\limits_{a}w_a
\psi_{\bar{\sigma}}(x+e_a \Delta t,t)e_a
\end{equation}

Where $\bar{\sigma}$ indicates the other fluid
component.$\psi_{\bar{\sigma}}$ and $\psi_{\sigma}$ are potential
functions and are taken as densities in this example.\\

A repulsive interactive force asks for a positive $G$ value. Large
$G$ value lead to sharp, non-diffusive interface, however give low
numerical stability at the same time.\\

Surface forces are incorporated into the multicomponent model in a
similar way as multi-phase model:

\begin{equation}
F_{surface}^{\sigma} (x,t)= -G_{surface}^{\sigma} \rho(x,t) \sum w_a
S(x+e_a \Delta t,t)e_a
\end{equation}




\section{Free Energy Model}
Swift et al (1995) developed a new thermodynamically consistent
binary fluid LB model by introducing an equilibrium state associated
with a free energy functional, corresponding pressure tensors and
chemical potentials. Correct choice of collision rules ensure that
the system evolutes toward to minimization of the free energy
functional.

The thermodynamic properties of a binary fluid system can be
described by a Laudau free energy functional.

\begin{equation}
\Psi=\int_V (\psi_b+\frac{\kappa}{2}(\partial_{\alpha} \phi)^2)dv
+\int_S \psi_s ds
\end{equation}

Where $\psi_b$ is the bulk free energy density and has a form:
\begin{equation}
\psi_b=\frac{c^2}{3} \rho ln\rho +A (-\frac{1}{2}\phi^2+
\frac{1}{4}\phi^4)
\end{equation}

$\phi$ is the order parameter represents the concentration of
components with the definition as:

\begin{equation}
\phi=\frac{\rho_a-\rho_b}{\rho_a+\rho_b}
\end{equation}


This free energy density function can lead to a phase separation
with $\phi=\pm 1$. The gradient term represents an energy
contribution and is related with surface tension. The second
integral is over the surface and describe the interactions between
fluids and the solid surfaces. $\psi_s$ is taken as:
\begin{equation}
\psi_s=-h \phi_s
\end{equation}

$\phi_s$ is the order parameter value at the solid surface.
Minimizing the free energy shows that the gradient in $\phi$ at the
solid surfaces is :

\begin{equation}
\kappa \partial_{\perp} \phi=\frac{d \psi_s}{d \phi_s}=-h
\label{wettingFE}
\end{equation}

The contact angle is related with the parameter $h$ and has a
formulation as:
\begin{equation}
h=\sqrt{2\kappa A} sign(\frac{\pi}{2}-\theta)
\sqrt{cos(\frac{\alpha}{3})[1-cos(\frac{\alpha}{3})]}
\end{equation}

Where $\alpha=cos^{-1}(sin^2 \theta)$.

Pressure tensor and chemical potentials are defined in an usual way:

\begin{eqnarray}
&\mu=\frac{\partial \Psi}{\partial \phi}=A(-\phi+\phi^3)-\kappa
\nabla^2 \phi&\\
&P_{\alpha
\beta}=(P_b-\frac{\kappa}{2}(\partial_{\gamma}\phi)^2-\kappa \phi
\partial_{\gamma \gamma}\phi)\delta_{\alpha \beta}+
\kappa(\partial_{\alpha}\phi)(\partial_{\beta}\phi)& \\
&P_b=\frac{c^2}{3}\rho+A(-\frac{1}{2}\phi^2+\frac{3}{4}\phi^4)&
\end{eqnarray}

The hydrodynamics and thermodynamics of the binary fluids are
described by the Navier-Stokes equations and a convection-diffusion
equation:

\begin{equation}
\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho u)=0
\end{equation}

\begin{equation}
\frac{\partial (\rho u)}{\partial t}+\nabla \cdot (\rho uu)=-\nabla
p_b + \nu \nabla ^2 (\rho u)+ \nabla[\lambda \nabla \cdot (\rho u)]
\end{equation}

\begin{equation}
\frac{\partial \phi}{\partial t}+\nabla (\phi u)= \mu \nabla ^2 \mu
\end{equation}

Where $\lambda= \nu (1-3\frac{c_s^2}{c^2})$.

A new distribution function $g_i(x,t)$ is introduced to describe the
concentration $\phi=\sum\limits_i g_i$ and is related to the
convection and diffusion. The distribution function $f_i(x,t)$ is
related to the fluid density and momentum as usual.

The time evolution equations using MRT scheme can be described as:
\begin{eqnarray}
\text{Collision Step:} & f'_i(x,t)=f_i(x,t)-M^{-1}SM(f_i-f^{eq}_i) & \\
& g'_i(x,t) =g _i^{eq}(x,t)& \\
\text{Streaming Step:} & f_i (x+e_i \Delta t,t+\Delta t) =f'_i(x,t)&
\\
&g_i (x+e_i \Delta t,t+\Delta t) =g'_i(x,t)&
\end{eqnarray}

Appropriate choice of equilibrium distribution fuctions can
reproduce the macroscopic equations in continum limit. The
equilibrium distribution functions should satisfy following
constrains:

\begin{eqnarray}
&\sum\limits_i f_i^{eq} = \rho& \\
&\sum\limits_i f_i^{eq}e_i = \rho u& \\
&\sum\limits_i f_i^{eq}e_i e_i=P + \rho uu & \\
&\sum\limits_i g_i^{eq} = \phi& \\
&\sum\limits_i g_i^{eq}e_i = \phi u& \\
&\sum\limits_i g_i^{eq}e_i e_i=\Gamma \mu I +\phi uu & \\
\end{eqnarray}

\subsection{Wetting boundary condition:}
Bounce Back rule is used to implement the non-slip condition. A
wetting boundary condition is given by equation (\ref{wettingFE})An
appropriate concentration value is assigned to solid nodes to
satisfy the equation(\ref{wettingFE}). $\nabla^2 \phi$ can be
calculated the same way as internal nodes. An example is give in
Figure () %FIGHER!!!!!!!!!!!!!!!!! and algorithm!!!!!!!!!

\section{Color Gradient Model}


An immiscible fluid model developed from Lattice Gas Cellular
Automata was introduced by Gunstensen and Rothman (1991) The
particles in this model are colored either red or blue, therefore it
is normally called color gradient method. The surface tension is
introduced by adding a perturbation to the collision operator while
keep the adherence to the Navier-Stokes equations in homogeous
regions. A recoloring step is involved after surface tension
perturbation calculation in order to achieve zero diffusivity of one
color into the other.

%====================================================================
A general color gradient lattice Boltzmann model can be summarized
as:



We use a D2Q9 lattice model as an example to illustrate the
algorithm of color gradient LBM. An extension to three dimensional
model is straightforward.

We use $f_i^r$, $f_i^b$ and $f_i$ denote the distribution functions
of a red fluid, a blue fluid and their combination respectively. The
densities and velocities can be calculated as:

\begin{eqnarray}
\rho^r=\sum\limits_i f_i^r & \rho^b=\sum\limits_i f_i^b &
\rho=\rho^r+\rho^b \\
&\rho u=\sum\limits_i e_i(f_i^r+f_i^b)&
\end{eqnarray}

An order parameter similar with the free energy model is defined:
\begin{equation}
\phi=\frac{\rho^r-\rho^b}{\rho^r+\rho^b}
\end{equation}

A regular MRT collision step is carried out for $f_i$:
\begin{equation}
f'(x,t)=f(x,t)-M^{-1}SM(f-f^{eq})
\end{equation}

Local kinematic viscosity $\nu$ is computed as a function of the
order parameter $\phi$
\begin{equation}
\nu=\nu_b+\frac{\phi+1}{2}(\nu_r-\nu_b)
\end{equation}

A perturbation is computed to generate the a surface tension.The
surface tension can be expressed as a local anisotropy in the
pressure: the pressure measured normal to the surface is larger than
that tangential to the surface (62) %%CITING!!!!!!!!!!!
The pressure in a LBM is proportional to the density, so surface
tension can be generated by placing preferentially particles in
directions normal to the interface rather than tangential. Mass and
momentum should be conserved.A color gradient G is defined as:

\begin{equation}
G(x,t)=\sum\limits_i e_i \sum\limits_j(f_j^r(x+e_i\Delta
t,t)-f_j^b(x+e_i\Delta t,t))
\end{equation}

Perturbation of the populations are:
\begin{equation}
f''_i(x,t) = f'_i(x,t) +A|G(x,t)|cos(2\theta)
\end{equation}

Where $\theta$ is the angle between the vector $e_i$ and color
gradient $G$. The parameter A is proportional to the magnitude of
surface tension $\sigma$ (D.Rothman and S.Zaleski Lattice Gas
Cellular Automata Cambridge University Press. Cambrige 1997):

\begin{equation}
\sigma=-\frac{192 \rho A}{\lambda}
\end{equation}

A redistribution of color enforce particles to move towards to the
regions occupied by particles with same color. This recoloring step
enable us to achieve the separation of the fluids. It is carried out
by the following maximization problem (Gunstensen and Rothman 1991)

\begin{equation}
W(f_i^{r''},f_i^{b''})=\max\limits_{f_i^{r''},f_i^{b''}}[\sum\limits_i
(f''_{r_i}-f''_{b_i})e_i]
\end{equation}

Subject to constrains:

\begin{eqnarray}
&\rho''_r=\sum\limits_i f_i^{r''}=\rho_r&\\
&f_i^{r''}+f_i^{b''}=f_i&
\end{eqnarray}

This maximization problem makes the color flux $H=\rho_{r}u_r-\rho_b
u_b$ have the same direction of color gradient and conserve the
local mass.

Many algorithms were developed to solve this maximization problem.
We use an algorithm proposed by T\"olke et al (2001). Comparing with
other existing algorithms, it does not generate any velocity
fluctuation for a plain interface and have a good numerical
stability. However the separation of phases is not as strong as the
other schemes.The algorithm is described as:

\begin{equation}
f_i^r=r_i^r+\frac{e_i \cdot G}{|G|} \ast min (f_i^b, f_{i+1}^b,
f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_{i+1}^r=r_{i+1}^r-\frac{e_i \cdot G}{|G|} \ast min (f_i^b,
f_{i+1}^b, f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_i^b=r_i^b-\frac{e_i \cdot G}{|G|} \ast min (f_i^b, f_{i+1}^b,
f_i^r, f_{i+1}^r)
\end{equation}

\begin{equation}
f_{i+1}^b=r_{i+1}^b+\frac{e_i \cdot G}{|G|} \ast min (f_i^b,
f_{i+1}^b, f_i^r, f_{i+1}^r)
\end{equation}

A general color gradient lattice Boltzmann model can be summarized
as:

\begin{itemize}
\item Single phase collision.
    %\begin{equation}
    %f_i=f_i^r+f_i^b
    %\end{equation}
\item Add a surface tension perturbation to $f'_i$ obtaining $f''_i$
\item Recoloring
\item Steaming
\end{itemize}



\section{Verification}
\subsection{Numerical results}
A series of drops surrounded by another immiscible fluid with same
viscosity are simulated to computed the surface tension. The
simulation results are similar as the SCMP results which confirm the
multi-component model agree with Laplace law.\\

Simulations of two fluid having different viscosities in a channel
is carried out and compared with theoretical predictions. An analytical
solution for a sharp interface is :

\begin{equation}
u_x(y)=\begin{cases}
\frac{1}{2\nu_1}G(L^2-a^2)+\frac{1}{2\nu_2}(a^2-y^2)& 0 \le |y| \le a\\
\frac{1}{2\nu_1} G(L^2-y^2) & a \le |y| \le L
\end{cases}
\end{equation}

Where $L$ is the width of the channel and $a$ is the location of
initial interface. The domain is periodic in x direction and bounded
by non-slip walls at $y=0$ and $y=128$ It is initialized with
substance 0 in the middle and substance 1 on both sides near the
walls. Initial density value 1 is applied to substance 0 and 1. 

The numerical results from Shan-Chen model is showed in Figure \ref{SC_v}: 

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{Shan-Chen_velocity.eps}
\end{center}
\caption{Simulated and analytical velocity profiles. The viscosity ratio between two components is 3, the density ratio is 1}\label{SC_v}
\end{figure}


The overall profile is as expected however it does not give good
agreement with the analytical solution. Density distributions are
plotted in Fig \ref{SC_d}

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{Shan-Chen_density.eps}
\end{center}
\caption{Density profiles of Component 1 and Component 0 with the initial density 1.0}\label{SC_v}
\end{figure}
A concentration drop in component 0 and a growth in component
1 are found. This phenomenon is due to the big diffusion on the
interfaces. The diffusivity can be reduced by increasing the
interaction forces however the numerical stability will be lowered
at the same time. Low diffusivity binary simulations are difficult
for Shan-Chen multi-component model due to this numerical stability
problem and deserve future exploration. It's worth mentioning that,
good agreements with analytical solution are obtained for same
viscosity immiscible fluids simulation. 

The same numerical tests are carried out using the Free Energy Model and Color Gradient Model. A viscosity ratio of 100 is introduced in the simulations. The densities of two components are still keep the same.   

\begin{figure}[h]
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{CG_poisseuille.eps}
%\caption{Simulated velocity of Color Gradient Model \label{CG_v}}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{FE_poisseuille.eps}
%\caption{Simulated velocity of Free Energy Model \label{FE_v}}
\end{minipage}
\caption{Simulated velocity of the Color Gradient Model (left) and the Free Energy Model (right) with viscosity ratio 100} \label{CGFE_v}
\end{figure}
	

The simulation results of both the Free Energy Model and the Color Gradient Model 
give excellent agreements with the analytical solution. The
concentration profile is also as expected.  Two immiscible fluid are
constrains in their respective areas and no unexpected diffusion is
discovered. We should notice that the interface thickness of the Free Energy Model is around $6$ lattice
units and the interface location is found at $x=38$ instead of their
initial position of $x=35$ A interface movement of 3 lattice is discovered from the result. The interface thickness of the Color Gradient Model is around 4, which is smaller than that of the Free Energy Model. A interface shifting is observed in color gradient model in
around 1.5 lattice units which is significantly smaller than the
free energy model.The recoloring algorithm separates two immiscible fluids very well
and gives nearly 0 diffusivity as we expected.

The maximum viscosity ratio for the Free Energy Model and the Color Gradient Model containing only static interfaces is around 120.




% \subsection{Capillary filling dynamics}
% Since Shan-Chen multi-component model cannot deal with binary fluids
% with different viscosities, the following comparisons are carried
% out between color gradient model and free energy model.
% 
% A capillary tube filled with binary fluids with viscosity contrast
% is simulated. A hydrophilic inner surface is associated with the
% tube. The energy liberated in wetting the surface is used to drive
% the fluid up the tube.
% 
% An analysis of capillary penetration is given by Pooley et al
% (2009). The length of the liquid column that penetrates the
% capillary could be calculated by surface tension $\sigma_{lg}$,
% capillary tube width $H$, dynamic contact angle $\theta$ and liquid
% viscosity $\eta$:
% 
% \begin{equation}
% l=(\frac{\sigma_{lg}H cos
% \theta}{3\eta})^{\frac{1}{2}}(t+t_0)^{\frac{1}{2}}
% \label{Washburnequation}
% \end{equation}
% 
% A system of $512 \times 32$ lattice units with periodic boundary
% condition in x direction is used in verification. Non-slip wetting
% boundary conditions are involved at the upper and lower sides of the
% tube.The simulation setup is showed in Fig( ) %FIGURE!!!!!!!!!!!
% An equilibrium gradient at the boundary $\partial_n
% \rho|_{bc}=-0.144$ is chosen for free energy model such that the
% equilibrium contact angle is $\theta=60^\circ$. The interface is
% initialized at $l=50$. Capillary filling distance are compared
% against with the Washburn equation (\ref{Washburnequation}).
% 
% %%FIGURE!!!!!!!!!!!!!
% 
% The distance of binary fluids along with capillary as a function of
% time. circle: $\nu_l=$, $\nu_g=$. square solid and dashed line are
% theoretical predictions.
% 
% We should notice that the maximum viscosity ratio of both free
% energy model and color gradient model decrease from around 120
% (static interface) to around 20 for dynamical interface simulation.




\subsection{Fingering}    
Capillary fingering is a kind of hydrodynamic instability which exists in various displacements during oil/gas production.
 One fluids is pushed by another one, with different viscosities, along a
channel with non-slip wall. A body force $G$ is applied on the
fluids. A growing finger of the driving fluid will be produced if the capillary
number $Ca$ is big enough.


\begin{equation}
Ca=\frac{u_t \rho \nu_2}{\sigma}
\end{equation}
Where $u_t$ is the velocity of the tip of the finger, $\nu_2$ is the
viscosity of driving fluid. $\sigma$ is the surface tension.

Fingering is investigated by two lattice Boltzmann binary models: the Free Energy Model and the Color Gradient Model.The standard bonce back boundary condition is imposed. Initially, the channel is occupied by component 1, component 2 is introduced at the begining of the channel, and is specified a constant velocity $u_t$. At the other end of the channel, a constant pressure boundary condition is imposed. A grid of (512 X 64) and a viscosity ratio of 10 are used in the simulations. Density of two components are still kept the same. The velocity at the tip of the fingers are measured by the code.




Figure () shows the evolutions of fingers simulated by the Free Energy Model. From up to the bottom, they are finger evolutions for surfance tension **** respectively. No finger will be produced if the surface tension is high. When the surface tension decrease, fingers are observed. However, the smaller surface tension, the less stable interfaces were produced. The forming of fingers is also related with the velocity at the tip of the fingers. Figure () shows the evlolutions of fingers simulated by the Color Gradient Model with the same surface tension but higher driving velocities. Fingers are formed in all simulations. 

Halpern and Gaver \cite{halpern_1994} studied the fingering phenomenon in a channel.  The width of the fingers produced by different capillary number $Ca$ are measured.Their results are ploted in the figure with a solid line(). A related width value with the width of the channel was used. Our results from the Free Energy Model and the Color Gradient Model are shown in the same figure. 

Good agreements are achived although some small discrepancies are found. These discrepancies might come from the boundary conditions.These numerical examples shows that, both the Free Energy Model and the Color Gradient Model are capable of simulating the capillary fingering phenomenon. It is worth mentioning that, the maximum viscosity ratio for dynamic interfaces simulations is around 20. Knowing this value can help people to identify if the Free Energy Model or the Color Gradient Model is suitable for their specific case.


\subsection{Snap-off Phenomena}

Snap-off is a type of mechanism for the fluid displacement in porous media.In a water wetting porous media which is comprised with pores and throat, the water will accumulate in some regions near the solid boundaries until the oil does not have contacting with the solid. Water then fill the throats and separate the oil into droplets. It takes place in the oil migration,and is of great industrial value. LBM has a strong capability of dealing with multi-component flows, therefore, I studied this important mechanism in pore scale fluid flow.


 \begin{figure}[!htbp]
\begin{center}
\includegraphics[width=3in]{snap-off.eps}
\end{center}
\caption{A pore/Throat geometry that used in the simulation}\label{snapoffgeo}
\end{figure}


A geometry showed in Figure(\ref{snapoffgeo}) is established to simulate the snap-off phenomena in pore scale. It is a tube comprises equilateral triangles with different incircle radius. The solid is set to water wetting and with a contact angle of $60^\circ$. A two dimensional view of the geometry is given in Figure (), due to the wetting property, the corner will be occupied by water. This is the region where water accumulate. A binary fluids system which comprises water of density and oil of density $1.0$ and viscosity of $0.1$ and $0.01$ repectively is simulated. I use a Free Energy Model with a grid 240x40x40, surface tension parameter $\kappa=0.04$, contact angle $20^\circ$ and a body force in X direction of $10^-5$. Initially, the water and oil occupy half of the tube respectively. The body force is applied to simulate the imbibition processing in the porous media.The oil will be displaced by the water.Figure() shows this displacement process. Water is founded to accumulate in the corners, and invaded the pores. 

As the oil was displaced by the water, the radii of curvature decreased. The capillary pressure determined by the radii of curvature decreased . As the radii of curvature decreases, the capillary pressure increases. The radii of the curvature kept decreasing while water was invading the throats. It broke the contacting of oil between pores if the surface tension was high enough. Some oil will stay in the pore, and will not be displaced by the water.  Figure () shows clearly this snap off phenomena. Water is represented by the red fluid, and oil is represented by the blue fluid. Distribution of oil and water is plotted in a interval of 1000 computational time steps. A part of the boundary is removed to better illustrate this processing. 



Surface tension is the critical condition for the emerge of snap-off. It should reach a critical value, then the oil will be cut into two part. From Laplace Law, we can find that, the surface tension is controlled by radii of curvature. It is determined by contact angle and geometry of the pore structure. A simulation with surface tension $\kappa=0.02$ was carried out to investigate the affection of surface tension on the snap-off.

 

We increase the minimum incircle radius of the geometry. This will decrease the maximum surface tension at which the oil phase can reach. The simulation results are showed in Figure (). We find that, no snap-off is produced.   

The contact angle is changed to $30^\circ$


\bibliographystyle{ieeetr}
\bibliography{ref}

\end{document}
